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Abstract. One of the most subtle points in the modern relativistic models for microarc- 
second astrometrical observations is the treatment of the influence of translational motion 



CN 
> 

of graving bo dl es on the light propagation. This paper descnbes numerical simulations 
of the light propagation in the gravitational field of moving gravitating bodies as well as 



summarizes the underlying theory. The simulations include high-precision numerical inte- 



grations of both post-Newtonian and post-Minkowskian differential equations of light prop- 
CD 

""■■■^ ■ agation and a detailed comparison of the results of the numerical integrations with various 

available approximate analytical formulas. The simulations has been performed both for hy- 
pothetical bodies with various parameters of trajectories as well as for all the major bodies 
of the solar system using the JPL ephemeris DE405/LE405 to calculate their motion. 



It is shown that for the accuracy of ~ 0.2 //as it is sufficient to use the well-known 
solution for the light propagation in the field of a motionless mass monopole and substi- 
tute in that solution the position of the body at the moment of closest approach between 



the actual trajectory of the body and the unperturbed light path (as it was first suggested 



by 



Hellings 



1986)). For a higher accuracy one should use either the post-Newtonian so- 
lution for uniformly moving bodies IQione^^JConeiki^^92 ) or the post-Minkowskian 



solution for arbitrarily moving bodies l Kooeikin &Schafer 1999). For astrometric obser- 



vations performed from within the solar system these two solutions guarantee the accuracy 
of ~ 0.002 yuas and are virtually indistinguishable from each other. 

Key words, astrometry - reference systems - relativity - gravitational lensing 

1. Introduction 



It is widely know n that extremely high accuracy of the future space-born astrometric mis- 
sions like GAIA ESaI bood IPerrvman et al. I boOlt Bienavme & Turon I 12002) and SIM 
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Shao 



1998) makes it necessary to formulate the reduction model of positional observations 



in a form fully consistent with General Relativity Theory (GRT). The relativ istic models o f 



posit ional observations has bee n form u lated by several groups of a uthors: 



Brumberg, Klioner, & 



de Felice et al. 



(2000), 



Kopeikin 



(1990) 



de Felice et al. 



Klioner & Kopeikin 



(2001), 



Klioner 



Klioner 



de Felice et al. 



1989), 



1998) 



J1992J) 

2003a). This paper is devoted to an 



investigation of one subtle point in any microarcsecond relativistic model of positional observa- 
tions. Namely, the influence of the translational motion of gravitating bodies on the light propa- 
gation is investigated here in great detail. 



After the pioneering work of 



Hellings 



( 1986) where it was suggested to compute the posi- 



tions of the gravitating bodies at the moment of closest approach between the body and the unper- 
turbed light ray and substitute these positions into the well-known solution for the light propaga- 



formulate more rigorous approaches to the problem 


Klioner Ill989 


1991 


Klioner & Kopeikin 


1992 


Kopeikin & Schafer 


1999; 


Kopeikin & Mashhoon 


2002 


). Detailed historical overviews 


can be found in Introduction of 


Kooeikin & Schafer 


Il999l and in Section 6 of 


Klioner 


(2003a! 



(see also 



Klioner 



l|2003bj». 

In this paper we perform extensive numerical simulations aimed at clarifying the ability of 
various approximate analytical formulas to reproduce the gravitational light deflection in the field 
of the solar system at the level of 0.1 - 1 //as as required by GAIA and SIM. 

Possible ways to compute the light trajectory in the gravitational field of moving bodies are 
summarized in Section |2] Section [5] explains general layout of the simulations. The results of 
the simulations are discussed in Section^] The suggestions for practical relativistic modeling of 
high-accuracy positional observations depending on the goal accuracy are formulated in Section 
13 In the Appendices we summarize the theoretical formulas concerning the influence of the 
translational motion of the gravitating bodies on light propagation. Most of these formulas are 
used in the simulations. Appendix 1X1 contains general formulas for null geodesies in the weak- 
field approximation. Both post-Newtonian and post-Minkowskian equations are given there. The 
most important theoretical results for light propagation in the post-Newtonian approximation are 
given in Appendix 151 The equations of light propagation in the post-Minkowskian approxima- 
tion are discussed in Appendix Icl The two point boundary problem for analytical solutions is 
discussed in AppendixlDl 



2. Possible ways to calculate the light propagation in the field of moving 
bodies 

According to the theory described in the Appendices there are several ways to calculate the light 
trajectory in the gravitational field of moving mass monopoles: 



1 . Numerical integration of the post-Minkowskian differential equations of light propagation 
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2. Analytical post-Minkowskian solution JC.28> -( IC.33> for arbitrarily moving bodies (the so- 
lution for the position of the photon JC.28I contains an integral that can be computed numer- 
ically, or estimated to be negligible for a particular purpose and thus neglected). 

3. Numerical integration of the post-Newtonian equations of light propagation dB.5> — (IB. 121 
with initial conditions JB . 1 31 — IB . 1 4i . 

4. Analytical post-Newtonian solution for uniformly moving mass monopoles (IB. 151 — (IB.26I 
with two free parameters xao and va to be related to the actual trajectory of the body. 

As discussed in the Appendices below the post-Minkowskian approximation scheme deals 
with expansions in powers of the gravitational constant G. Velocities of gravitating bodies are 
not considered as small in the post-Minkowskian approximation scheme which is sometimes 
called weak-field fast-motion approximation. The first post-Minkowskian approximation implies 
that all terms of order 0(G 2 ) are neglected. The post-Newtonian approximation scheme operates 
with expansions in powers of c _1 . In the post-Newtonian approximation scheme velocities of 
gravitating bodies are considered to be small. This approximation scheme is sometimes called 
weak-field slow-motion approximation. In the first post-Newtonian approximation terms of order 
0(c~ 4 ) are neglected in the equations of light propagation. One can prove that in the case of light 
propagation the formulas of the first post-Newtonian approximation are linear in G and, therefore, 
contained in those of the first post-Minkowskian approximation (see e.g. AppendixO. 

The most accurate way to calculate the light propagation is clearly the first one, i.e. numerical 
integration of the post-Minkowskian differential equations of motion for the photon. The errors 
of those differential equations of motion come from the effects of the second post-Minkowskian 
(or post-Newtonian) approximation which are known to be negligible for the solar system appli- 
cations. Numerical integration of the post-Newtonian equations of motion for the photon dB.5i — 
JB. 12l > can also be used to calculate the light trajectory in the field of arbitrarily moving bodies. 
However, since the post-Newtonian equations of motion are contained in the post-Minkowskian 
equations of motion (AppendixO. the former can be used for an internal consistency check of 
the whole calculation rather than for an independent computation of the light path. Our numerical 
experiments show that for solar system applications the results of numerical integrations of the 
post-Newtonian and post-Minkowskian equations of motion are identical within the errors of the 
first post-Minkowskian approximation. 

The analytical post-Minkowskian solution given in Appendix lC.4l is surely the most accurate 
analytical solution for the problem. However, (1) the solution for the photon's position involves 
an integral which should be in principle computed numerically, and (2) the post-Minkowskian 
solution is relatively expensive as far as the computing time is concerned since it contains the 
retarded moment of time to be computed by iterations (see below). On the other hand, the full 
accuracy of the analytical post-Minkowskian solution is not necessary to attain the accuracy of 1 
yuas for the solar system applications. Simpler analytical solutions of sufficient accuracy can be 
found instead. 
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The fully analytical post-Newtonian solution given in Appendix IB .41 describes the light tra- 
jectory in the field of uniformly moving gravitating bodies having the coordinates 

XA(t)=XAO+VA(t-t ), (1) 

where xao and Va are two arbitrary constant vectors. These constants can be related to the actual 
trajectory of the body JC^ ph (f) in different ways with the hope that the errors related to the non- 
uniformity of the body's trajectory will be minimized in some sense. The principal goal of this 
paper is to check if the analytical post-Newtonian solution with some reasonable choice of the 
constants Xaq and Va can describe the gravitational light deflection with an accuracy of 0.1 - 1 
fias. 

For the analytical post-Newtonian solution one has to choose either a fixed point or a straight 
line as the model trajectory Xa(0 of the body which should be distinguished from the actual 
trajectory of the body x e ^ h (t). In this paper we consider six choices for the constants xao and 
Va in the post-Newtonian solution: four solutions for a body at rest named Pi, Pi, P3 and P' 3 
and two solutions for a body moving with constant velocity L\ and L2 (see Figure 0. Each of 
the considered solutions uses the actual position (and possibly velocity) of the body at some 
reference moment of time t = f re f . The most simple solution Pi uses the fixed position jc^ ph (f ) 
of the body at the moment of observation f re f = t . The fixed position *^ ph (f ca ) of the body at the 
moment f re f = f ca of closest approach between the body and the unperturbed light ray is used for 
the second solution Pi. The moment f' d can be calculated as 



: max 



t e ,t - max 



g=V ~ - x e * b (t ), 
c 



( g-(x p (t )-xf(t )) 
c\g\ 2 



(2) 
(3) 



where t e is the moment of emission of the photon, and x p (t) is the (unperturbed) trajectory of the 
photon. If the source is situated outside of the solar system one can put t e = -00 and the outer 
max can be omitted in @. The retarded moment of time t* is used as f re f in the solution P3. This 
moment of time is defined by 



1 



t* + -\x p (t)-x A p \t*)\ = t. 



(4) 



The moment f is relatively expensive to calculate since the equation @ is an implicit one and 
one has to use some kind of iterations to solve it (e.g., Newton's method). That is why, one can 
try to substitute t* by its simplified version f which can be directly calculated 



f =t--\x p (t)-x t f(i)\. 



(5) 



The solution P' 3 uses the position :c^ ph (?*') of the body at f re f = ?*'. 

The solutions L\ and L2 are the solutions for a body moving with a constant velocity. In these 



two cases the parameters xao and va of ( IB. 151 — JB.26I are chosen to coincide with the actual 
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solution 


tat 




v A 


Pi 


to 







Pi 




xT<.n 





P 3 


f 


*?(<*) 





P' 


t" 


xf(f) 






u t xTito) tT«o) 

Table 1. Choice of the constants for the six analytical post-Newtonian solutions for the light 
trajectory considered in the present paper (see text for further explanations). 



light path 




Fig. 1. Actual trajectory x e ^ (?) of body A and the five simplified trajectories used to model the 
light propagation within the post-Newtonian approximation scheme: three fixed positions P\ 
(xA(t) = x e ^ h (t ) = const), P2 (xA(t) = *^ ph (? ca ) = const) and P3 (*a(0 = ^ ph (f*) = const), and 
two trajectories with constant velocity L\ (jca(0 = x^ ph (f fl ) + i^ ph (f c )) (f - t a ) and L2 (xaO) = 
jc^ ph (f ca ) + i^ ph (f ca )) (f - f ca ). The position P', mentioned in the text is not shown on the sketch. It 
is situated on the actual trajectory of the body close to P3. 

position and velocity of the body at the moments t m t = t and f re f = f ca , respectively (therefore, 
the trajectories of the bodies used in L\ and L2 are tangents to the actual trajectory of the body at 
these two moments of time). The choice of the trajectories XAif) are summarized in Tableland 
Figure^ 

3. Simulations 

It is clear that numerical integration of the differential equations of light propagation can be 
performed only for sources situated at some finite distance from the gravitating body (the end 
point of numerical integrations is anyway at some finite distance since it is defined by the position 
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light path 




n 



body's trajectory 




Fig. 2. Four vectors appearing in the calculations: vector fx is the unit light direction at the point 
of emission x p (to), n is the unit light direction at the point of observation x p (f), k is the unit 
coordinate direction from x p (to) to x p (t), and cr is the unit direction of the light propagation for 
t — > —oo. Formal definitions of these vectors are given in AppendixiDl 



of observer). Light propagation from a source at a finite distance to the observer represents a 



two point 



Klioner 



boundary problem for the differential equations of light propagation. As discussed by 



(2003a) the goal of the relativistic reduction of observations in this case is to relate the 



unit direction n of the light propagation at the moment of observation to the unit direction k from 
the point of light emission x p (to) to the point of light observation x p (t) (see, Figure|2j. 

In each individual simulation we fix the points of emission and observation for all the solu- 
tions. These points are computed by numerical integration of the post-Minkowskian differential 
equations, so that for this most accurate solution we compute vectors n, fx, x p (to), x p {t) and k 
with the maximal possible accuracy. Using the formulas given in Appendix iDl we then solve 
the two point boundary problem for all the analytical solutions discussed in our simulations and 
compare the vector n from the post-Minkowskian numerical integration to the vectors n from the 
analytical solutions. We choose the distance between the point of emission and the gravitating 
body sufficiently large so that the differences between the vectors n do not change (within an 
accuracy of 0.001 /vas) when we further increase the distance. Therefore, we can claim that our 
results are valid also for sources at infinity. Moreover, the simulations have shown that the dif- 
ferences between vectors n for the points of emission at lower distances are smaller than those 
for sources at infinity. This means that the maximal errors given below are valid also for sources 
situated at finite distances from the observer (e.g. for the solar system objects). 

In order to test the analytical formulas in different situations and check internal consistency 
of the simulations we have done three independent series of simulations with different choices 
for the trajectories x e ^ h (t) of the bodies. For all the simulations the masses of the gravitating 
bodies have been taken from the JPL ephemeris DE405/LE405 while the radii and other parame- 



ters of the bodies were taken from Weissman. McFadden & Johnson 



( 1999). The three series of 



simulations can be described as follows. 
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A. Parabolic trajectories with constant acceleration (Table|5Jl. For each trajectory the velocity 
and the acceleration of the body at the moment of closest approach between the body and the 
photon coincide with the maximal possible barycentric velocity and acceleration of the corre- 
sponding body of the solar system. The impact parameter for each trajectory is chosen so that at 
the moment of closest approach the distance between the photon and the body is equal to the ra- 
dius of the corresponding body (expect for the three lines in Table|2]with fixed minimal allowed 
angular distance if/mm between the gravitating body and the direction of light propagation as seen 
by the observer; the minimal avoidance angle for each of these three bodies (Earth, Sun and 
Moon) are calculated from the condition that the minimal Sun avoidance angle is iffi. = 35°). 
The distance between the observer (satellite) and the gravitating body is taken to be the maximal 
possible distance between the GAIA satellite and the corresponding body of the solar system 
(some of the effects under study become larger with increasing this distance and it is important 
to use the maximal possible distance for the simulations). For the calculation of that maximal 
distance the satellite was supposed to be situated exactly at the Lagrange point L2 of the Earth- 
Sun system. After fixing all these parameters the initial conditions of the photon trajectory are 
still not unique and this freedom can be used to check all possible mutual orientations of the 
velocity vector of the photon and the velocity and acceleration vectors of the body at the mo- 
ment of closest approach. The directions of all these three vectors are independent of each other. 
Routinely, in a coordinate system where the direction of the velocity of the photon is fixed we 
check 50 uniformly distributed directions for each of the two other vectors to find the maximal 
value of the effects under investigation. In several cases we have checked that a finer grid of 
mutual orientations does not lead to any changes in the maximal differences between the models 
given in Table|2] 

B. Circular coplanar trajectories (Table |3J. The observer is supposed to be situated exactly 
at the Lagrange point L2 of the Earth-Sun system. Each gravitating body moves along a circular 
orbit with realistic semi-major axis and mean motion. All the orbits are coplanar. All possible 
configurations of the observer and the gravitating body have been checked with a step of 0.01 
of the siderial period of the corresponding body. The impact parameter of the light ray (i.e. the 
minimal distance between the photon and the body) is chosen to be not smaller than the radius 
of the body, but can be larger to meet the requirement imposed by the minimal Sun avoidance 
angle ilf. = 35°. For all the bodies for which the minimal Sun avoidance angle influences their 

fc ~ mm & 

observability (Sun, Mercury, Venus, Earth and Moon), the simulations has been done with and 
without taking the minimal Sun avoidance angle of i^® in = 35° into account. For each mutual 
configuration of the body and the observer 36 different initial positions and velocities of the 
photon have been chosen so that in the view plane of the observer the observed directions of 
light are uniformly distributed around the observed position of the gravitating body. Table [3] 
contains maximal differences between the models for all the light trajectories investigated for 
each of the body. 
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C. Realistic trajectories on the basis of the JPL ephemeris DE405 (Table 0}. Simulations 
similar to B have been performed using the JPL numerical ephemeris DE405/LE405 for the 
trajectories of the gravitating bodies. The orbit of the observer is taken to be a realistic Lissajous- 
type orbit ab out the Lq. point of the Earth-Sun system and was computed using the algorithm 



suggested by 



Mignard 



(2002). For each body the minimal Sun avoidance angle of 35° has been 
taken into account while choosing the parameters of the light rays. All mutual configurations of 
the observer and the gravitating bodies between 1 January 2008 and 3 1 December 2020 have been 
checked with a step of 1 day (the total time span covers 4749 days). The impact parameters of the 
light rays are taken in the same way as for simulation B. Also in the same way as in simulation 
B we have checked 36 different directions of the light ray uniformly distributed relative to the 
line connecting the observer and the gravitating body. We also checked in several cases that 
increasing the number of the observed light directions from 36 to, say, 360, does not change the 
maximal differences between the models given in Table[3] 

For simulations A, B and C an ANSI C program has been written. Since the effects we are 
looking for can be as small as 10~ 3 ;uasx: 5 • 10" 15 , it is not sufficient to perform the computations 
using standard "double" 64-bit arithmetic which provides only 16 decimal digits. Routinely, we 
have used 80-bit arithmetic on an Intel processor (18 decimal digits). Some additional accuracy 
checks have been performed using 128-bit arithmetic on an Ultra-Sparc processor (34 decimal 
digits) and have shown that the 80-bit arithmetic does not produce substantial roundoff errors and 
is sufficient for our purposes. Note that since the simulations took several weeks of computing 
time on a Pentium III processor running at 600 MHz (about one million photon trajectories for 
each gravitating body were checked), it was hardly possible to perform the simulations within 
a reasonable time using a software environment emulating arithmetic with arbitrary precision 
(Maple, Mathematica, etc.). 

For the numerical integrations of both the post-Minkowskian and post-Newtonian differen 
tial eq uations of motion we have used the Everhart integrator described, e.g. in 



Everhart 



(1974 



1985). In our program the Everhard integrators of the orders 7, 11, 13, 15, 19, 23 and 27 are 
implemented, and the internal coefficients of the integrator are coded with an accuracy consis- 
tent with the 128-bit arithmetic. This makes it possible to perform the numerical integration with 
very high precision (at least, up to 34 decimal digits) in a quite efficient way. Our investigation 
showed that the number of internal iterations within the integrator can be chosen to be as low as 
1 (or at most 2) without any loss of the resulting accuracy. This can be understood as a conse- 
quence of almost straight trajectories of the light. The integrator of order 19 was found to be the 
most efficient for our calculations. The final global accuracy of the numerical integrations was 
controlled by integrating the solution backwards. The integration is repeated automatically with 
higher accuracy parameters of the integrator until the required goal accuracy is reached. 

The most time-consuming part of the calculations is the numerical integration of the post- 



Minkowskian differential equations of light propagation given in Appendix IC.2I . In order to 
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maximally speed up the numerical integration the code calculating the right-hand side of JC.11> 



has been optimized using the Maple package CODEGEN dChar et al. I I 1 993 >. This optimization 



reduced the number of necessary float-point operations roughly by a factor of 2. 



4. Results of the simulations 

The maximal differences between the vectors n from the numerical post-Minkowskian solution 
and the various analytical solutions are summarized in Tables The following conclusions 
from the numerical simulations can be formulated: 

1 . The results of the three independent simulations are in good agreement with each other. This 
serves as an internal consistency check of the simulations. 

2. The post-Minkowskian analytical model is virtually indistinguishable from the numerical 
integration of the post-Minkowskian differential equations of light propagation and leads to 
errors of order 0.002 //as. 

3. As expected the naive model for the light propagation involving the post-Newtonian analyti- 
cal solution for the body being at rest at its position at the moment of observation (the model 
Pi) is too inaccurate and leads to errors exceeding 1 mas. Note that our software code for 
the model Pi does not check if the formally calculated impact parameter of the light rays 
exceeds the radius of the body (such a check is, of course, easy to implement, but if the 
software detects that the formal impact parameter is smaller than the radius of the body, this 
can serve only as an evidence of the inaccuracy of the model P\ and the situation cannot be 
corrected within the model Pi). This is the reason why the errors for the model P\ given in 
Tables fJJ0] sometimes exceed the maximal possible gravitational light deflection due to the 
corresponding body. 

4. The post-Newtonian analytical solution for the body being at rest at its position at the moment 
of closest approach (P2) or at the retarded moment of time (P3) are virtually indistinguishable 
from each other for the solar system applications (e.g., for Jupiter the maximal difference of 
these models does not exceed 7.5 x 10~ 4 //as). 

5. Any of these two models (P2 and P3) allows one to attain an accuracy of ~ 0.18 //as for the 
realistic trajectories of the gravitating bodies (Table[4}. The errors of ~ 0.75 //as appearing in 
the simulation with parabolic trajectori es also follow f rom a simplified analytical considera- 



tions and are quoted e.g. in Table 1 of Klioner 



(2003a) as upper estimates of the effects. The 



reason for the discrepancy between these estimates fo r the parabolic tr ajectories and those 



Klioner 



(2003a, p. 1590, above Eq. 



for the realistic trajectories were already discussed by 
(34)) where the realistic values for Jupiter and Saturn were predicted. One can see that the 
predicted realistic values and the values from Tableware in good agreement. 
6. Simplified calculation of the retarded moment of time as given by (model P' 3 ) increases 
the errors in the light deflection. The errors attains ~ 0.3 juas for the realistic trajectories of the 
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gravitating bodies. Additional simulations show that the following formula gives the retarded 
moment with sufficient accuracy (this represents one Newtonian iteration for Eq. 0) 



'"'='- , , ep h , A ' P =X p (t)-xT(t). (6) 

c \P\ ~ *A W ' P 

A solution P'l which is similar to Pi , but with f le f = f" has the same errors (within the level 
of 0.001 fias) as the solution P3 with exact value of the retarded moment f, e f = f . 

7. The post-Newtonian analytical model for a body moving with a constant velocity is indistin- 
guishable from the post-Minkowskian model within the accuracy of 0.002 juas provided that 
the position and velocity of the body on the rectilinear model trajectory of the body coincide 
with the actual positions and velocities of the body at the moment of the closest approach 
(model L2). 

8. If the position and velocity of the body on the rectilinear trajectory coincide with the actual 
positions and velocities at the moment of observation (model L\ ) the error exceeds 0. 1 



Let us also note that increasing the minimal Sun avoidance angle i^® in would reduce the errors 
for the Sun, Mercury, Venus, Earth and the Moon (all the gravitating bodies which are closer to 
the Sun than the observer). Therefore, the numbers in Tables |2j0] can be considered as upper 
estimates for if/^ > 35°. 



5. Concluding remarks 



The resu l ts of ou r numerical simulations are in good agreement with the theoretical discussion by 



Klioner 



( 2003a). These results allows one to formulate the following practical recommendations 



for data processing of microarcsecond positional observations performed from the solar system: 



- If an accuracy of 0.2 fias is sufficient one can employ the simple post-Newtonian analytical 
model for the light propagation in the gravitational field of a motionless body and substitute in 
that model the actual position of the body evaluated either at the moment of closest approach 
f ca or at the retarded moment f . If f is used, Eq. (|6j can be employed to calculate it with 
sufficient accuracy. 

- If an accuracy below 0.2 //as is required one should use either the full post-Minkowskian 
analytical model as given in Appendix IC.4I (neglecting the integral g(t^, t) in dC.30» or the 
post-Newtonian analytical model for a body moving with a constant velocity (Appendix lB.4> 
and choose the constants Xao and Va of that model to coincide with the actual position and 
velocity of the body at f ca . 



body 



S.A. Klioner, M. Peip: Numerical Simulations of the Light Propagation 1 1 

6 S 2 pM pN,Pi pN,P 2 pN,P 3 pN,P' 3 pN,L\ pN,L 2 



Sun* 


1.76 • 10 6 


15.0 


17.7 


43.3 


17.8 


17.8 


17.8 


17.7 


17.7 


Sun, i^ =35° 

1 ~ min 


13600 


0.001 


0.001 


0.003 


0.002 


0.002 


0.002 


0.001 


0.001 


Mercury* 


82.9 


0.0 


0.0 


93.2 


0.016 


0.016 


0.154 


0.601 


0.0 


Venus 


493 


0.0 


0.0 


812 


0.058 


0.058 


0.178 


0.357 


0.0 


Earth* 


574 


0.0 


0.0 


8.08 


0.058 


0.058 


0.058 


0.0 


0.0 


Earth, i// mil> = 35° 


3.64 


0.0 


0.0 


0.001 


0.001 


0.001 


0.001 


0.0 


0.0 


Moon* 


25.9 


0.0 


0.0 


1.44 


0.003 


0.003 


0.003 


0.0 


0.0 


Moon, if/ mir[ = 19° 


0.086 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


Mars 


116 


0.0 


0.0 


143 


0.010 


0.010 


0.058 


0.096 


0.0 


Jupiter 


16300 


0.001 


0.002 


26900 


0.746 


0.746 


0.847 


0.292 


0.002 


Saturn 


5780 


0.0 


0.0 


86100 


0.196 


0.196 


0.250 


0.036 


0.0 


Uranus 


2530 


0.0 


0.0 


5960 


0.047 


0.047 


0.110 


0.085 


0.0 


Neptune 


2080 


0.0 


0.0 


6530 


0.050 


0.050 


0.106 


0.081 


0.0 



Table 2. The result of simulations using parabolic trajectories x e ^ (f). All quantities are given in 
yuas. The second column 6 is the maximal light deflection angle for the body computed with the 
post-Minkowskian differential equations of light propagation, while 6 2 is the approximate value 
of the post-post-Minkowskian terms which are neglected in our considerations. The maximal 
errors of the post-Minkowskian analytical solution (pM) and the six post-Newtonian solutions 
(pN) are given. The meaning of P, and L, is discussed in Section [5] and summarized in Table 
^ The number "0.0" means that the actual value of the effect is less than 0.001 pas. The bodies 
designated with a star, e.g. "Mercury*", cannot be observed by GAIA because of the minimal Sun 
avoidance angle of at least 35°. The minimal impact angle tf/ m \ n for the Sun and the Earth coincide 
with the minimal Sun avoidance angle of t^° in = 35°, while that for the Moon is approximately 
calculated from simplified orbits of GAIA and the Moon. 



Appendix A: Equations of null geodesies 

Here we summarize the formulas for the null geodesies in a weak gravitational field used in this 
paper. The metric tensor g a p is supposed to have the form 

gap = r\afs + Kp, (A.l) 
where 

i7 a p = diag(-l,+l,+l,+l) (A.2) 
is the flat Minkowski metric and h a p = h a p{t, x) is the non-Galilean part of the metric. 
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S S 1 pM pN,Pi pN,P 2 pN,P 3 pN,P' 3 pN,Li pN,L 2 



Sun* 


1.75 • 10 6 


14.8 


17.5 


33.3 


17.5 


17.5 


17.5 


17.5 


17.5 


Sun, if/ s " n = 35° 


12900 


0.001 


0.001 


0.002 


0.001 


0.001 


0.001 


0.001 


0.001 


Mercury* 


83.1 


0.0 


0.0 


315 


0.013 


0.013 


0.091 


0.192 


0.0 


Mercury, iffi* = 35° 


0.007 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


Venus* 


493 


0.0 


0.0 


33000 


0.058 


0.058 


0.156 


0.156 


0.0 


Venus, if, Sun = 35° 

' min 


493 


0.0 


0.0 


33300 


0.058 


0.058 


0.145 


0.144 


0.0 


Earth* 


574 


0.0 


0.0 


13.8 


0.0 


0.0 


0.0 


0.0 


0.0 


Earth, ^ Slln = 35° 

' ' min 


3.85 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


Moon* 


26.2 


0.0 


0.0 


3.06 


0.001 


0.001 


0.001 


0.0 


0.0 


Moon, </f Slm = 35° 

' r mm 


0.092 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


Mars 


116 


0.0 


0.0 


257 


0.006 


0.006 


0.033 


0.022 


0.0 


Jupiter 


16300 


0.001 


0.002 


21300 


0.139 


0.139 


0.202 


0.035 


0.002 


Saturn 


5780 


0.0 


0.0 


31300 


0.020 


0.020 


0.035 


0.008 


0.0 


Uranus 


2080 


0.0 


0.0 


3550 


0.003 


0.003 


0.009 


0.003 


0.0 


Neptune 


2530 


0.0 


0.0 


3700 


0.002 


0.002 


0.007 


0.003 


0.0 



Table 3. The maximal difference between the models in the simulations with circular coplanar 
trajectories jc^ p (f) for the gravitating bodies. All quantities are given in /ias. See caption of Table 
|2and Section|3]for further details. 



A.1. Post-Newtonian approximation 

First, we use the standard post-Newtonian assumptions on the orders of magnitude of h a p with 
respect to the formal parameter c _1 : 



hoo = 0(c- 2 ), 
h 0i = 0(c- 3 ), 

hj = 0(c~ 2 ). (A.3) 
Then the equations of null geodesies in the first post-Newtonian approximation read 



x' = I c 2 h m - hoot x k x' - (hit i - ^ Aw ) x k x' -\ h m x' 



2 *w,( ■ "Wit, " i ■ -i is., i 2 a-m i 2 
- 1^ hokj ~ h jk )j x J x k x' -c {h m - h k,i) ** 

-h iKt ^+0{c- 2 ), (A.4) 

where for any small latin index i except for t one has Aj — A and A_, = j f A. Note that the 
velocity of the photon x k = O(c). Eq. \AA\ is valid in the post-Newtonian approximation which 



body 
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Sun 


13000 


0.001 


0.001 


0.002 


0.001 


0.001 


0.001 


0.001 


0.001 


Mercury 


0.012 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


Venus 


493 


0.0 


0.0 


30000 


0.058 


0.058 


0.147 


0.149 


0.0 


Earth 


4.75 


0.0 


0.0 


0.001 


0.0 


0.0 


0.0 


0.0 


0.0 


Moon 


0.125 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


Mars 


116 


0.0 


0.0 


385 


0.008 


0.008 


0.040 


0.025 


0.0 


Jupiter 


16300 


0.001 


0.002 


19600 


0.175 


0.175 


0.255 


0.038 


0.002 


Saturn 


5780 


0.0 


0.0 


27900 


0.030 


0.030 


0.052 


0.008 


0.0 


Uranus 


2080 


0.0 


0.0 


3550 


0.004 


0.004 


0.014 


0.003 


0.0 


Neptune 


2530 


0.0 


0.0 


3700 


0.002 


0.002 


0.008 


0.003 


0.0 



Table 4. The maximal difference between the models in the simulations with realistic trajectories 
x^ ph (f) of the gravitating bodies taken from the JPL planetary ephemeris DE405/LE405. All 
quantities are given in //as. See caption of Tableland Section[3]for further details. 



is also called weak-field-slow-motion approximation. The assumption /zo, = 0(c~ 3 ) means that 
the velocity of the gravitating bodies are considered to be small with respect to c. Therefore, 
these equations cannot be used if the velocities of the gravitating bodies are large. 

A. 2. Post-Minkowskian approximation 

Second, we use the post-Minkowskian properties of the metric tensor and consider all the com- 
ponents of h a p to be of first order with respect to gravitational constant G: 



h 00 = 0(G), 
h 0i = 0(G), 

hj = 0(G). (A.5) 
No expansion in terms of c _1 is used here. Then, one has 



x l = ^ c 2 /ioo,/ - hoo,k x* x l - yi ikJ - i hujj 

(1 1 \ , k , 

- 1 - hokj - h jkA x J x x -c (h i,k - hokj) x' 



h iKt x k -ch m ,, + 0(G 2 ). (A.6) 



The only formal difference between jA.6> and jA.4> is the last term -c /io,,/ which has first order 
with respect to G, but is 0(c~ 2 ) and has been omitted in ( IA.4t . Eq. ( IA.6> has been derived without 
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any assumption on the velocity of the bodies. This equation is valid in the first post-Minkowskian 
approximation which is also called weak-field limit. 

A.3. Initial-value problem for the equations of motion 

Initial value problem for the differential equations iAA\ and JA.6t can be formulated as 



X(t ) = Xq, 

x(t ) = cns(t ), fi-fi = l. (A.7) 
Function s can be generally determined from the condition of the geodesic to be a null one: 

gaP iehP = Q, ** = (i,i±j, (a.8) 



(A.9) 



which for metric ( IA. li gives 

_ V(l - /zoo) (1 + hjjjj'/j.j) - hom' 
l+hijfi'fii 

Note that s is a function of time and position (via the metric components h a p) as well as of the 
direction fi, but for a given trajectory of the photon it can be considered as a function of time. 

Taking into account the orders of magnitude of h a p given by iA.3\ in the first post-Newtonian 
approximation one gets 



* = 1 - \ hoo - l - hij/j^ - hoin 1 + 0(c- 4 ). (A.10) 

It is easy to see that this expression is valid also in the first post-Minkowskian approximation 
(that is, it contains all terms of dA.9l > linear with respect to G). 



Appendix B: Light propagation in the post-Newtonian approximation 

B.1. Metric tensor 

The non-Galilean components of the metric tensor in the post-Newtonian approximation read 



2 

%)= w(t,x) + 0(c' 4 ), 
4 

hoi= — ^w'(t,x) + 0(c- 5 ), 

h u = -| Sij wit , x) + 0(c- 4 ) . (B . 1 ) 

For moving mass monopoles the potentials have the form 

V GMa ru ™ 

A A 



z 



GMa 
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where r A = \r A \, 
r A (t,x) = x-x A (t), 
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(B.3) 

(B.4) 



x A (t) is the position of body A, and x A (f) is the velocity of body A. 

The metric (IB. 1 1 with potentials ( IB.2I — (IB.3i coincides w ith the metric tensor in the 



IAU 



(2001) (see also 



Soffel at al. 



(2003) 



Barycentric Celestial Reference System given by the 
for a discussion). The higher -order multipole moments of the bodies' gravitational fields ( IB. 2 
JB.3i also discussed in liAirlikooil) will be ignored in this paper. 



6.2. Differential equations of light propagation 

The post-Newtonian equations of motion of a photon then read 



-Z 



GM A 



(A A n A +B A v + C A v A ) + 0(c 2 ), 



(B.5) 



A A = 2 + y — 45, 

B A = 4(l ~ a)S - (1 -/3)(2 + y), 

C A = -4(l-a), 



(B.6) 
(B.7) 
(B.8) 



a=\-n A v, (B.9) 

p = l-n A -v A , (B.10) 

y=l-v-v, (B.ll) 

8=1 -vv A , (B.12) 



where n A = r A /r A , v = jc/c, Va = JCa/c 
6.3. Initial-value problem 

From the general formulas of Section lA~3l the initial-value problem for ( IB.5> - (IB.12> reads 

X(t ) = Xq, 

x(t ) = cns(to), HH=l, (B.13) 
and 



s(t) = 1 - 




(B.14) 
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B.4. Analytical solution for a body in uniform motion 

Two analytical solutions of Eqs. (IB.5I — (IB. 121 and (IB.13I — JB. 141 are known: 1) the classical 
solution for a body at rest, i.e. Xa = const; and 2) a solution for a body moving with a constant 
velocity. The first solution is clearly contained in the second one. An approximate analytical 
solution for bodies having a constant velocity, t hat is for x A (t) = Xao + Va (t - to), where Xaq 



const and Va - const, has been first derived by 



Klioner 



(1989) (see also 



Klioner & Kopeikin 



1992). The solution reads 



x(t) = x(t ) + cfi s(t ) (t - t ) + Ax(t , t) 

-Ax(to) (t - 1 ), 

- x(t) = fi s(to) + - Ax(t) - - Ax(to), 
c c 



c 

Ax(t , t) 

8a 
Ia 



- Ax(t) ■ 

c 



Ia 
Ja 



■ fi x (r A o x g A ), 
■fi-v A , 



1 



1 



\g A W r A\- gA' r A \g A W r A()\- gA' r A0 

\gA\\rA\ + g A -r A 



Ja = log 

\g A W r Ao\ + g A -rA0 

r A (t) = x(t) - x A (t), 



r A o = x(t ) - x A (to), 



■z 



2GM A 



\g A \ 



d A -i A +g A -j A ]+0(c-*), 

c c 



\rA\(\g A \\ r A\- g A - 'a) 
\gA\ 

ir*r 



(B.15) 
(B.16) 

(B.17) 
(B.18) 
(B.19) 
(B.20) 

(B.21) 
(B.22) 
(B.23) 
(B.24) 

(B.25) 
(B.26) 



This solution has two parameters Xao and Va - Note that the values of Xao an d v a are arbitrary 
and it is apriori unclear how to relate Xao to the actual trajectory JC^ ph (f) of the body in order to 
minimize the errors (see Sectional. 

The classical solution of ( IB . 5I >— dB~T2l > for a body at rest can be easily restored from ( IB. 151 — 



JB .26l > by setting Va = that implies g A - H and \g h \ = 1. 



Appendix C: Theoretical results in the post-Minkowskian approximation 

C.1. Metric tensor 



In the post-Minkowskian approximation the metric t ensor o f a system of a rbitrarily 



monopoles can be written (see 



Kopeikin & Schafer 



(1999); 



moving mass 



Kopeikin & Mashhoon 



(2002)) with 



retarded potentials simila r to the L ienard-Wiechert potentials which are well known from the 



classical electrodynamics (Jackson I1974I) . The metric reads 
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h m = \<W(t,x) + 0{G 2 ), 
h 0i = ~'W i (t,x) + O(G 2 ), 

h tj = ^W ij (t,x) + 0(G 2 ), (Cl) 



Tvc^Y^^-r; 1 ), (c.2) 
«u"fex) = y^r.cv;'. (c.3) 
tv%, x) = Y ^ r; 1 + 2 r. v A ; v A j ) , (c.4) 

where r * = * - x* A , r* = |r*|, x* A = x A (<*), v* = £x A /c, n A = r*/r*, & = 1 - n* • v*, 

/ -.-1/2 _ 

r» = ( 1 - • v^J . Here and below in Section|CJthe position Xa, velocity Xa and acceleration 
Xa of the gravitating bodies are computed at the time moment t* A = f A (t, x) being the retarded 
moment of time defined by the following implicit equation 

f*+ir*=*. (C.5) 

Note that r A (t) = r A J3 t + 0(r 2 A lc 2 ), so that formally TV = w + 0(c" 2 ), TV' = tv' + 0(c- 2 ), 
<W ij = £' 7 w + <9(c~ 2 ), and Eqs. (IB. 1MB. 4 can be easily restored from dC.U - dC.4t within the 
first post-Newtonian approximation. 



C.2. Differential equations of light propagation 
Substituting dC.H into dA.6> one gets 



x = TV,- - 4 TV* x k x l - 4 TV* ,x k x l + 4 TV*' ,i* i' 

- 4 TV , i* + 4 TV * , x j x k x l + 4 TV * . i J ' i* i'' 

c ' c c 

+4w vi*-4'**'*.-**-4'w*f** 

c 2 ' c 2 • c 2 

+4 r W , ', + <5(G 2 ). (C.6) 

c z 

Here, as usual, A ; , = Jj A is the partial derivative with respect to t, and for any latin index i 
with except for t one has A >( - = 4^ A. Since the potentials TV, TV' and TV'-' are given by Eqs. 
dC.2> - dC.4> explicitly as functions of the retarded time f A - f A (t, x) one has to use 



d dA(t*x) dt* 

_ A( „ )= _J__d, (C.7) 

»^_£^ (C8 , 

ox' dt. ox' ox 1 
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where from dC.5> 



df 

-4=IK\ (C9) 
at 

g = - c -V»r- (cio) 



Substituting (IC.2I — dC74l > into \C.6l and using dC.7t -( lC8t one gets the explicit form of the 
equations of motion of a photon in the first post-Minkowskian approximation 



GMi r 



*2 r/3 



I- 



- (&a n A +S A v + C A v A + D A a* A ) 



+0(G 2 ), 



(CM) 



M A = [r; 2 y-26l] rf (r; 2 + e«) 



(C12) 



—2 



+ 26* \/3 t 6* - e, (2a, - S t )} + 4f./J„(a, - <?,) 
+77,/?, (r; 2 y - 25*(2a* - 5,)) , 
C A =r; 4 [4tf.or.-£.y] 



+2r; 



0. {2e,a, -p t 6,}-2£,p*a. 



+4 77, a»j8, 6„, 
D A =4Tf a t p, 6 t c~ l r* A , 



(C.13) 



(C14) 
(C.15) 



a t = 1 - n* A ■ v, (C.16) 

p. = l-nl-v* A , (C.17) 

y=l-v-v, (C.18) 

6» = l-v-v\, (C.19) 

£.=(«:-i;)c" 1 4 (c.20) 

f^Cai-v)^ 1 ^, (C.21) 

r}, = (a A -v^c-'rl, (C.22) 



and v = c- 1 | x, a* = ^v* = c" 1 ^ *^). 

Note that if the accelerations a* A are considered to be of gravitational nature (and, therefore, 
a* A - 0(G)) the terms in ( IC.1U containing a* A are of order 0(G 2 ) and can be discarded in the first 
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post-Minkowskian approximation. However, in our theory no equations of motion of the gravi- 
tating bodies are considered, so that the accelerations a* A can be of any nature and are considered 
to be given. That is why we retain the accelerations a* A in jC.ll> in the first post-Minkowskian 
approximation. 

From Eqs. dC.lH - dC.22> one can restore the post-Newtonian equations of motion dB.5i — 
JB. 12l > neglecting in ( IC. 1 II — dC.22l all the terms formally of order 0{c~ 2 ) (e.g., all the 
acceleration-dependent terms) and taking into account that in ( IC. 1 li — dC.22i the position of the 
bodies are calculated at the retarded moment of time f* while in dB.5> -( IB"T2l at moment t . The 
latter circumstance gives the following relations: r^(f) = r* A -r* A v* A +0(c~ 2 ), r A (t) = r* A /3*+0{c~ 2 ), 
n(t) = (n* - v*)/#, + 0{c 2 ). 

Since from dC.l 1> it follows that 

x = xo + cn(t-t Q ) + 0(G), ftfi=l (C.23) 

in the first post-Minkowski approximation in the right-hand side of Eq. dC.ll> one can put 

x = xq + cfi(t - t ), fi fi-1 (C.24) 

and, therefore, v — fi, which implies, e.g., y = 0. Formally, this makes Eq. JC. 1 li integrable in 
quadratures (since the right-hand side is simply a function of time). We prefer here not to do so 
and retain Eq. dC.ll> in the form of a differential equation as given by dC.6> . 

C.3. Initial-value problem 
Initial value problem for Eq. dC.ll> 



X(t ) = Xq, 

x(t ) = cfil(t ), nn= 1 
with 



A "A 



(C.25) 

(C.26) 
(C.27) 



C.4. Analytical solution in the post-Minkowskian approximation 

Th e explicit form of the analy tical solution in the first post-Minkowskian approximation derived 



by 



Kopeikin & Schafer 



(1999) reads 



x(t) = x(t ) + c/i s(f ) (t - t ) + Ax(t , t) 
-Ax(to) (t - to), 



(C.28) 
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- x(t) = n s(to) + - Ax(t) - - Ax(to), 
c c c 



Ax(t , t) = - J] [f(t) - f(t ) + g(t , t)j + 0(G 2 ), 

I a x (r* x u) 
fit) = tAb. / - (/i - vi) log(rl or.) 



««b,o= r 4V (itf,-e.>7,)/ix(jiix/i) 

+a» log(r* a.) (77, (/i - v* A ) - V~ 2 c~ l r* A a* A ) 



cdt, 



iA±(i) = -2 



2GM A T, 0» / n A 1 x ("1 x 

c 2 



+(2-0 t )^-2v;J + O(G 2 ). 



(C.29) 
(C.30) 
(C.31) 



(C.32) 



(C.33) 



Here g(to, t) are integrals depending on the accelerations of the bodies. The integrals can be 



transformed into integrals with respect to retarded time f A as described in 



Kopeikin & Schafer 



( 1999). In the right-hand sides of dC.30t and dC.33t one should put x(t) = x(to) + c/x(t- to). 

Note that the solution for - Ax(t) is exact within the first post-Minkowskian approximation 
(i.e. all the omitted terms are of order G 2 ). It is remarkable that this solution can be written in 
such a compact closed form for arbitrary motion of the gravitating bodies. Note also that the 
solution depends only on the position and velocity of the bodies in a single moment of time - the 
retarded moment t A corresponding to time t - and does not depend on the previous motion of the 
bodies. If all the terms which are formally of order 0(c~ 4 ) (these are the terms at least quadratic 



with respect to v^) are dropped both in dC.33t and JB.24t the two formulas are equivalent. One 
can also show that Eqs. jC.30t - dC.3U agree with Eqs. ( IB. 17> - dB.2H in the first post-Newtonian 
approximation provided that v A is considered to be constant. 

Taking into account that for any function A' and x(t) — Xq + c fi(t - to) 



1 n o 

- A\t, x(t)) = - A% x(t)) + — A\t, x(t)) cn } 
at at ox 1 

= ? ^zA\f,x{t)) + ^- AXf,x{t)) Ci i', 
/J* at ox 1 



(C.34) 



one can prove that a time derivative of dC.33> exactly coincides with the post-Minkowskian 
equations of motion for a photon ( IC. 1 li — ( IC22I provided that v is taken to be equal to fi in the 
latter (this is allowed within the first post-Minkowskian approximation as discussed above). 



A different way to derive the Kopeikin-Scha fer solution 



moving with constant velocities is given by 



Klioner 



C"2"51 - dc331 for the case of bodies 



(2003b) who has shown how to derive this 



solution combining the post-Newtonian solution (IB. 151 — dB.26i for a motionless body (v A = 0) 
with a suitable Lorentz transformation. 
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Appendix D: Two point boundary value problem for the analytical solutions 

Let us consider analytical equations of light propagation in the form valid for both analytical 
models discussed above 



x(t) = x(t ) + cfi s(t ) (t - t Q ) + Ax(t , t) 

-Ax(t )(t-t ) + O(e 2 ), (D.l) 
1 11, 

- x(t) = fi s(t ) + - Ax(t) - - Ax(t ) + 0(e 2 ), (D.2) 
c c c 

with 



Ax(t , to) = 0, 
lim - Ax(t) = 0. 

t— » — CO Q 



(D.3) 
(D.4) 



These two conditions can be proven to be valid for both analytical models considered above. The 
small parameter e can be identified with G for the post-Minkowskian solution and with c~ 2 for 
the post-Newtonian solution, so that in both solutions Ax(to, t) = 0(e) and c~ l Ax(t) = 0(e). 
Let us define several vectors 



(T = lim - x(t), 



n = n(t) 



x(t) 

W)\' 

R(t , t) 



* = * M= |*(^)f 
R =R(t ,t)=x(t) -x(t ). 



(D.5) 
(D.6) 

(D.7) 
(D.8) 



These definitions are illustrated on Figure [2] Vector er is the unit direction of the light path at 
past null infinity (t — > -oo). Vector n(t) is the unit direction of light propagation at time moment 
t (so that n(to) = ft and lim n(t) = <t). Vector k is the unit vector connecting the points x(ty) 

t—*—oo 

and x(t). Now, if one considers the points jc(fo) and x(t) as given constants, one gets from JD. II — 
JD.4> the following formulas approximately solving the two point boundary value problem for 
the differential equations of motion 



1 , 
o-=(i-fix(- Ax(t Q ) Xfi) + 0(e L ), 
c 

s(t) = 1 + - fi ■ Ax(t) + 0(e 2 ), 



n — fi + fi x 

k — fi + fix 

+0(e 2 ), 



- Ax(t) - - Ax(t ) 
c c 



— Ax(to) + — Ax(t ,t) 
c \R\ 



xfi) + 0(e 2 ), 



(D.9) 
(D.10) 
(D.ll) 

(D.12) 
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n = k + k x 



- Ax(t) - — Ax(t , t) 
c \R\ 



x k 



+0(e 1 ). (D.13) 

Note that jD.9> follows from ( ID.1H for t — > -oo, that JD. 1 31 is a combination of ( ID.1U and 
( ID.12L and that ( ID.10> agrees with the expressions IB. 14b and dC.26> for s(t) given above for 
the post-Newtonian and the post-Minkowskian solutions, respectively. From (ID. 121 and the cor- 
responding expressions for x(t) and Ax(to, t) one can see that for a given impact parameter of the 
light ray the difference between vectors k and fi becomes smaller for greater distances between 
the gravitating body and point of emission x(to). 

Note that Ax(to, t), Ax{t) and Ax(to) on the right-hand sides of dD.9> -( ID.13> depend on fi. 
Formally, considering only analytical orders of magnitude one could replace fi by k in these 
formulas. However, this works well only when the impact parameters computing for the unper- 
turbed trajectories with directions fi and k are sufficiently close to one another. This is not always 
the case. For example, this assumption is wrong in the gravitational lens limit. It is more reliable 



to calculate vector fi from given k by a numerical inversion of (ID. 12b and then calculate n from 
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